// eee.Sheffield.PZ.Imaging
//
// Copyright ?Ping Zou, 2007
// sg71.cherub@gmail.com

namespace eee.Sheffield.PZ.Imaging
{
    using System;
    using System.IO;
    using System.Drawing;
    using System.Drawing.Imaging;
    using eee.Sheffield.PZ.Math;
    using eee.Sheffield.PZ.Imaging.ImageFilter;
    using System.Collections;
    using System.Collections.Generic;

    /// <summary>
    /// bifurcation search
    /// </summary>
    public class LiuSunBifurcationSearch
    {
        #region Fields
        private LiuSunTrackResult _trackResult;
        private List<PZPoint> _cpCandidateArray;
        private List<PZDirection> _vdCandidateArray;
        private List<double> _initrCandidateArray;
        private Bitmap _bifurcationImage;
        // parameters
        private double _Noff = 6.0;
        private double _Nskip;
        #endregion

        #region Properties
        public List<PZPoint> CentrePointCandidateArray
        { get { return _cpCandidateArray; } }
        public List<PZDirection> VesselDirectionCandidateArray
        { get { return _vdCandidateArray; } }
        public List<double> InitialRadiusCandidateArray
        { get { return _initrCandidateArray; } } 
        #endregion

        #region Constructor
        public LiuSunBifurcationSearch(LiuSunTrackResult trackResult)
        {
            _trackResult = trackResult;
            _cpCandidateArray = new List<PZPoint>(0);
            _vdCandidateArray = new List<PZDirection>(0);
            _initrCandidateArray = new List<double>(0);
        } 
        #endregion

        public void Apply(Bitmap srcImage)
        {
            _bifurcationImage = new Bitmap(srcImage);
            // search along left edge
            SearchAlongLeftEdge(srcImage, _trackResult, ref _cpCandidateArray, ref _vdCandidateArray, ref _initrCandidateArray, ref _bifurcationImage);
            // search along right edge
            SearchAlongRightEdge(srcImage, _trackResult, ref _cpCandidateArray, ref _vdCandidateArray, ref _initrCandidateArray, ref _bifurcationImage);
            // save bifurcation image
            _bifurcationImage.Save("bifurcation.bmp", ImageFormat.Bmp);
        } // Apply()

        private void SearchAlongLeftEdge(Bitmap srcImage, LiuSunTrackResult trackResult, 
            ref List<PZPoint> cpCandidateArray, ref List<PZDirection> vdCandidateArray, 
            ref List<double> initrCandidateArray, ref Bitmap bifurcationImage)
        {
            int width = srcImage.Width;
            int height = srcImage.Height;
            int trackResultLength = _trackResult.Count;

            // search along left edge
            // 1. average gray-level value for pixels along left edge
            double GL = 0.0;
            for (int i = 0; i < trackResultLength; i++)
            {
                int x = (int)_trackResult.EdgePoint1Array[i].x;
                int y = (int)_trackResult.EdgePoint1Array[i].y;
                if (x >= 0 && x < width && y >= 0 && y < height)
                    GL += (double)srcImage.GetPixel(x, y).R;
            }
            GL /= (double)trackResultLength;

            // 2. for each centre point
            int trackResultLengthM1 = trackResultLength - 1;
            for (int i = 0; i < trackResultLengthM1; i++)
            {
                PZPoint pCi = _trackResult.CentrePointArray[i];
                PZPoint pLi = _trackResult.EdgePoint1Array[i];
                PZPoint pCiA1 = _trackResult.CentrePointArray[i + 1];
                PZPoint pLiA1 = _trackResult.EdgePoint1Array[i + 1];
                PZDirection di = _trackResult.VesselDirectionArray[i];
                PZDirection dLi = new PZDirection(di);      // normal to di, point to Left 
                dLi.Rotate(90);

                int distancepLiToiA1 = (int)(System.Math.Floor(pLi.Distance(pLiA1)));                
                for (int j = 0; j < distancepLiToiA1; j ++)
                    // for each edge point between pL(i) and pL(i + 1)
                {
                    // pL(i, j)
                    PZPoint pLij = new PZPoint();
                    pLij.ExtendFrom(pLi, di, j);
                    // pL+(i, j)
                    PZPoint pLijP = new PZPoint();
                    pLijP.ExtendFrom(pLij, dLi, _Noff);
                    // extract scanline through pL+(i, j)
                    int widthOfScanline = (int)(System.Math.Round(_trackResult.RadiusArray[i] * 2 + 4));
                    //int widthOfScanline = (int)(System.Math.Round(_trackResult.RadiusArray[i] * 3));
                    ExtractScanline exsl = new ExtractScanline(pLijP, di, 1, widthOfScanline);
                    exsl.Apply(srcImage);
                    PZMath_vector scanlineProfile = exsl.ScanlineProfile;
                    // detect branch
                    double G = 0.0;
                    for (int k = 0; k < widthOfScanline; k++)
                        G += scanlineProfile[k];
                    G /= (double)widthOfScanline;
                    if (G > GL)
                        // find a branch
                    {
                        // triangular matched filter to detect cp candidate
                        double initr = _trackResult.RadiusArray[i];
                        int halfWidth = widthOfScanline / 2;
                        SunMatchedFilter matchedFilter = new SunMatchedFilter();
                        PZMath_vector output = matchedFilter.TriangularMatchedFilter(scanlineProfile, (int)initr);
                        
                        // centre point candidate
                        int maxPos = output.MaxIndex() - halfWidth;
                        PZPoint cpCandidate = new PZPoint();
                        cpCandidate.ExtendFrom(pLijP, di, maxPos - halfWidth);

                        // vessel direction candidate, pC(i, j) -> cpCandidate
                        PZPoint pCij = new PZPoint();
                        pCij.ExtendFrom(pCi, di, j);
                        PZDirection vdCandidate = new PZDirection();
                        vdCandidate.GetDirectionFrom2Points(pCij, cpCandidate);

                        // radius candidate, use a Sun rectangular matched filter
                        PZMath_vector rectOutput = matchedFilter.RectangularMatchedFilter(scanlineProfile, (int)initr);
                        int centrePosition = rectOutput.MaxIndex() - halfWidth;
                        if (centrePosition <= 0)
                            continue;
                        double signalLevel = matchedFilter.SignalLevel(scanlineProfile, centrePosition, initr);
                        double backgroundLevel = matchedFilter.BackgroundLevel(scanlineProfile, centrePosition, initr);
                        double rollOffPoint = matchedFilter.RollOffPoint(signalLevel, backgroundLevel);
                        PZMath_vector edgePoints = matchedFilter.EdgePoints(scanlineProfile, rollOffPoint, ProfileShape.Peak);
                        double initrCandidate = (edgePoints[1] - edgePoints[0] + 1.0) / 2.0;
                        
                        // add candidates to the array
                        cpCandidateArray.Add(cpCandidate);
                        vdCandidateArray.Add(vdCandidate);
                        initrCandidateArray.Add(initrCandidate);

                        // plot on bifurcation image
                        bifurcationImage.SetPixel((int)cpCandidate.x, (int)cpCandidate.y, Color.Red);

                        // skip Nskip
                        _Nskip = halfWidth - 1;
                        j += (int)_Nskip;
                    }
                } // for j
            } // for i
        } // SearchAlongLeftEdge()

        private void SearchAlongRightEdge(Bitmap srcImage, LiuSunTrackResult trackResult,
            ref List<PZPoint> cpCandidateArray, ref List<PZDirection> vdCandidateArray, 
            ref List<double> initrCandidateArray, ref Bitmap bifurcationImage)
        {
            int width = srcImage.Width;
            int height = srcImage.Height;
            int trackResultLength = _trackResult.Count;

            // search along left edge
            // 1. average gray-level value for pixels along left edge
            double GR = 0.0;
            for (int i = 0; i < trackResultLength; i++)
            {
                int x = (int)_trackResult.EdgePoint2Array[i].x;
                int y = (int)_trackResult.EdgePoint2Array[i].y;
                if (x >= 0 && x < width && y >= 0 && y < height)
                    GR += (double)srcImage.GetPixel(x, y).R;
            }
            GR /= (double)trackResultLength;

            // 2. for each centre point
            int trackResultLengthM1 = trackResultLength - 1;
            for (int i = 0; i < trackResultLengthM1; i++)
            {
                PZPoint pCi = _trackResult.CentrePointArray[i];
                PZPoint pRi = _trackResult.EdgePoint2Array[i];
                PZPoint pCiA1 = _trackResult.CentrePointArray[i + 1];
                PZPoint pRiA1 = _trackResult.EdgePoint2Array[i + 1];
                PZDirection di = _trackResult.VesselDirectionArray[i];
                PZDirection dRi = new PZDirection(di);      // normal to di, point to right 
                dRi.Rotate(-1.0 * 90);

                int distancepRiToiA1 = (int)(System.Math.Floor(pRi.Distance(pRiA1)));                
                for (int j = 0; j < distancepRiToiA1; j ++)
                    // for each edge point between pR(i) and pR(i + 1)
                {
                    // pR(i, j)
                    PZPoint pRij = new PZPoint();
                    pRij.ExtendFrom(pRi, di, j);
                    // pR+(i, j)
                    PZPoint pRijP = new PZPoint();
                    pRijP.ExtendFrom(pRij, dRi, _Noff);
                    // extract scanline through pR+(i, j)
                    int widthOfScanline = (int)(System.Math.Round(_trackResult.RadiusArray[i] * 2 + 4));
                    //int widthOfScanline = (int)(System.Math.Round(_trackResult.RadiusArray[i] * 3.0));
                    ExtractScanline exsl = new ExtractScanline(pRijP, di, 1, widthOfScanline);
                    exsl.Apply(srcImage);
                    PZMath_vector scanlineProfile = exsl.ScanlineProfile;
                    // detect branch
                    double G = 0.0;
                    for (int k = 0; k < widthOfScanline; k++)
                        G += scanlineProfile[k];
                    G /= (double)widthOfScanline;
                    if (G > GR)
                        // find a branch
                    {
                        // triangular matched filter to detect cp candidate
                        double initr = _trackResult.RadiusArray[i];
                        int halfWidth = widthOfScanline / 2;
                        SunMatchedFilter matchedFilter = new SunMatchedFilter();
                        PZMath_vector output = matchedFilter.TriangularMatchedFilter(scanlineProfile, (int)initr);
                        
                        // centre point candidate
                        int maxPos = output.MaxIndex() - halfWidth;
                        PZPoint cpCandidate = new PZPoint();
                        cpCandidate.ExtendFrom(pRijP, di, maxPos - halfWidth);

                        // vessel direction candidate, pC(i, j) -> cpCandidate
                        PZPoint pCij = new PZPoint();
                        pCij.ExtendFrom(pCi, di, j);
                        PZDirection vdCandidate = new PZDirection();
                        vdCandidate.GetDirectionFrom2Points(pCij, cpCandidate);

                        // radius candidate, use a Sun rectangular matched filter
                        PZMath_vector rectOutput = matchedFilter.RectangularMatchedFilter(scanlineProfile, (int)initr);
                        int centrePosition = rectOutput.MaxIndex() - halfWidth;
                        if (centrePosition <= 0)
                            continue;
                        double signalLevel = matchedFilter.SignalLevel(scanlineProfile, centrePosition, initr);
                        double backgroundLevel = matchedFilter.BackgroundLevel(scanlineProfile, centrePosition, initr);
                        double rollOffPoint = matchedFilter.RollOffPoint(signalLevel, backgroundLevel);
                        PZMath_vector edgePoints = matchedFilter.EdgePoints(scanlineProfile, rollOffPoint, ProfileShape.Peak);
                        double initrCandidate = (edgePoints[1] - edgePoints[0] + 1.0) / 2.0;
                        
                        // add candidates to the array
                        cpCandidateArray.Add(cpCandidate);
                        vdCandidateArray.Add(vdCandidate);
                        initrCandidateArray.Add(initrCandidate);

                        // plot on bifurcation image
                        bifurcationImage.SetPixel((int)cpCandidate.x, (int)cpCandidate.y, Color.Red);

                        // skip Nskip
                        _Nskip = widthOfScanline - 1;
                        j += (int)_Nskip;
                    }
                } // for j
            } // for i
        } // SearchAlongRightEdge()
    }
}
